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Abstract 

The maximum-final-mass trajectory of a proposed 
configuration of the Advanced Launch System is pre- 
sented. A model for the two-stage rocket is given; the 
optimal control problem is formulated as a parame- 
ter optimization problem; and the optimal trajectory 
is computed using a nonlinear programming code called 
VF02AD. Numerical results are presented for the con- 
trols (angle of attack and velocity roll angle) and the 
states. After the initial rotation, the angle of attack goes 
to a positive value to keep the trajectory as high as pos- 
sible, returns to near zero to pass through the transonic 
regime and satisfy the dynamic pressure constraint, re- 
turns to a psotive value to keep the trajectory high and 
to take advantage of minimum drag at positive angle of 
attack due to aerodynamic shading of the booster, and 
then rolls off to negative values to satisfy the constraints. 
Because the engines cannot be throttled, the maximum 
dynamic pressure occurs at a single point; there is no 
maximum dynamic pressure subarc. 

To test approximations for obtaining analytical solu- 
tions for guidance, two additional optimal trajectories 
are computed: one using untrimmed aerodynamics and 
one using no atmospheric effects except for the dynamic 
pressure constraint. It is concluded that untrimmed 
aerodynamics has a negligible effect on the optimal tra- 
jectory and that approximate optimal controls should 
be able to be obtained by treating atmospheric effects 
as perturbations. 

List of Symbols 

Cl lift coefficient 
Cd drag coefficient 
C m pitching moment coefficient 
D aerodynamic drag (lb) 
g local gravitational acceleration (ft /sec 2 ) 
h altitude (ft) 
i orbit inclination (rad) 

I sp vacuum specific impulse (sec) 

J performance index 
l aerodynamic reference length (ft) 
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L aerodynamic lift (lb) 
m vehicle mass (slugs) 

M* aerodynamic pitching moment (ft lb) 

Af Mach number 
P penalty function 
p atmospheric pressure (lb/ft 2 ) 
q dynamic pressure (lb/ft 2 ) 

Sb aerodynamic reference area (ft 2 ) 

T thrust (lb) 

T vae vacuum thrust (lb) 
f, staging time (sec) 

V velocity (ft/sec) 

a angle of attack (rad) 

7 flight path angle (rad) 

6 thrust gimbal angle (rad) 

9 pitch angle (rad) 

A longitude (rad) 
p velocity roll angle (rad) 
p atmospheric density (slug/ft 3 ) 
r latitude (rad) 
r normalized time 

u rotational velocity of earth (rad /sec) 
ip heading angle (rad) 

Subscripts 

6 body axes 
eg center of gravity 
e exit 
/ final 
i inertial 
o initial 
8 sea-level 
tv wind axes 

I. Introduction 

A program is under way to develope an unmanned, 
all-weather, launch system for placing medium to large 
payloads (~ 120,0001b) into low-earth orbit. A prospec- 
tive design for this Advanced Launch System \Aw/ is 
shown in Fig. 1 to be composed of a core vehicle and 
a booster. Both the booster and the core are ignited 
at launch, and staging occurs when all the booster pro- 
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pellant is consumed. Payload mass can be increased by 
adding another booster. 

Part of the design process is to iterate the vehicle de- 
sign and trajectory design until a reasonable combina- 
tion is achieved. This paper is concerned solely with the 
optimal trajectory design of the proposed configuration. 
The objective is to find the trajectory which maximizes 
the final mass (since the engines burn throughout the 
trajectory, this is also a minimum final time problem). 
Any remaining propellant can be considered for conver- 
sion to payload or a decrease in launch weight. The 
physical model is that of flight over a rotating, spher- 
ical earth with an exponential atmosphere. Launched 
vertically from the surface of the earth, the payload is 
to be placed into perigee of an 80nm by 150nm transfer 
orbit. Because of structural considerations, there is a 
limit on the amount of dynamic pressure the vehicle can 
withstand. 

This study has had several goals: (a) to determine the 
maximum-final-mass trajectory of the proposed ALS, 

(b) to generate initial Lagrange multipliers for a shooting 
code to investigate neighboring extremal guidance, and 

(c) to determine if atmospheric effects (pressure thrust 
and aerodynamics) can be considered as a perturbation 
to vacuum thrust and gravity for guidance law devel- 
opment. While only (a) and (c) are reported here, (b) 
requires the use of an exponential atmosphere. Hence, 
the dynamic pressure limit based on a standard atmo- 
sphere has been lowered to have the same effect in an 
exponential atmosphere. 

In Section 2, a model is presented for the proposed 
ALS configuration. Then, the optimal control problem 
is formulated in Section 3 and converted into a param- 
eter optimization problem in Section 4. This is done 
for relative ease in obtaining an optimal trajectory. Nu- 
merical results are presented in Section 5 in the form 
of optimal controls, states, and dynamic pressure. Also 
contained in Section 5 are two additional optimal trajec- 
tories based on untrimmed aerodynamics and neglected 
atmospheric effects. Finally, conclusions are presented 
in Section 6 . 


II. Physical Model 

In this section, a physical model for the Advanced 
Launch System (ALS) is defined. It includes the equa- 
tions of motion for flight over a rotating, spherical earth 
with an exponential atmosphere and the mass, propul- 
sion, and aerodynamic properties of the vehicle. 

Equations of Motion 

Since sideslip causes drag, the vehicle is assumed to fly 
at zero sideslip angle, so that only the angle of attack 
gives the orientation of the vehicle relative to the free 
stream. The direction of the lift vector is then controlled 


through the bank angle or, more specifically, through the 
velocity roll angle. 

A three-degree-of-freedom model for vehicle motion 
can be obtained from a six-degree-of- freedom model by 
one of two aerodynamic approximations: untrimmed 
aerodynamics or trimmed aerodynamics. For a rocket, 
untrimmed aerodynamics is equivalent to setting the 
thrust gimbal angle to zero and ignoring the aerody- 
namic pitching moment. On the other hand, with 
trimmed aerodynamics, it is assumed that the pitch rate 
is zero (pitching moment equals zero) so that the gim- 
bal angle can be determined as a function of the angle 
of attack. 

In view of the above comments, the three-degree-of- 
freedom equations of motion relative to the earth are 
given by (Ref. 3) 
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In these equations, A is the longitude, t is the latitude, h 
is the altitude above mean sea level, V is the velocity, 7 
is the flight path angle, tp is the heading angle, m is the 
mass, r = r, + h is the distance from the center of the 
earth to the vehicle center of gravity, w is the angular 
velocity of the eath, D is the drag, L is the lift, T is the 
thrust, I , p is the specific impulse, 6 is the gimbal angle 
of the thrust vector, a is the angle of attack, and fi is the 
velocity roll angle. With regard to signs, a positive roll 
angle generates a negative heading toward the south. 

For trimmed aerodynamics, the pitching moment, 
which is the sum of the aerodynamic pitching moment 
and the thrust pitching moment, is set equal to zero, and 
the resulting expression solved for the thrust gimbal an- 
gle. With reference to Fig. 1 and by assuming that 6 is 
small, this process leads to 
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Figure 1: Force and Moment Nomenclature 


where Ma is the aerodynamic pitching moment ant l t 
is the distance from the center of gravity to the exit 
plane of the engines. Because 6 is dependent on the 
aerodynamic pitching moment and the moment is de- 
pendent on the pitching moment coefficient, it results 
that 6 is linear in a with the coefficients varying with 
time. Aerodynamics is discussed in further detail later 
in this section. 

Eqs. (1) have two singularities: V = 0 in the 7 and 
the \p equations and 7 = j in the ip equation. To re- 
move the V singularity and to clear the launch tower, 
the vehicle is flown vertically for 3 sec with the angle of 
attack and the bank angle being chosen so that 7 = 0 
and ip = 0. To remove the 7 singularity, the vehicle is 
pitched over at constant heading (ip = 0) for 1.0 sec at 
a constant negative pitch rate $ whose optimal value is 
determined. Since 6 = 7 + a, the angle of attack during 
pitch-over is given by 

a = J - 7 + *(* - 3) . (3) 

m 

Finally, the bank angle is chosen to make tp = 0. With 
a flat earth model, /i = 0. 

Earth 

The earth is taken to be a rotating, spherical body 
whose surface is described by the mean searlevel radius 
r, and whose gravitational acceleration varies with alti- 
tude according to the inverse-square law 

9 = 9.(7)’ <<> 

where g,r] represents the earth’s gravitational parame- 
ter. Sea-level gravitational acceleration g, } r,, and the 
rotational velocity of earth ui are known constants given 

33 ft 

r, = 2.09256725E+7 ft , g, = 32.174 ^ 

u = 7.2921 158E-5 — . (5) 


The atmosphere is represented by the exponential 
functions 


— =ez P {-A) . £ = ex P(j:'> 

p t Ai p, 


where the scale-height constants are given by 

Ai = 23,800 ft , A 2 = 23,200 ft 


( 6 ) 

(7) 


and the searlevel values of the density and pressure are 
p, = .002377 ^ , p, = 2, 116.24 ^ . (8) 

Finally, the speed of sound is given by 



where 7 = 1.4 is the ratio of specific heats of air. 


Mass Characteristics 

The ALS configuration consists of a core vehicle as 
depicted in Fig.l. The take-off mass of the ALS con- 
sists of the inert vehicle mass, the propellant mass, pay- 
load mass, payload margin mass, and the payload fairing 
mass (Table 1). 


Table 1: Mass Characteristics 


Vehicle 

Vehicle Component 

Take-off Mass 
(slugs) 

Core 

Inert Mass 
Propellant 
Payload 

Payload Margin 
Payload Fairing 

5,474.29 

45,974.38 

3,729.71 

372.97 

1,215.89 


■Ki&issaH 


Booster 

Inert Mass 
Propellant 





Core + 
Booster 

Toted Take-off Mass 

108,574.93 


The center of gravity is located relative to a coordinate 
system whose origin is at the tip of the core vehicle, 
whose x axis is down the longitudinal axis, and whose 
y axis is toward the booster. For the first stage, the 
vehicle center of gravity is assumed to have coordinates 

x c , = 165.45 ft , y cg = 10.36 - .0388/ ft (10) 
so that l t is constant and has the value 

110.81 ft (11) 


sec 
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where / = 276.26 ft is the length of the core vehicle. 
Actually, x cg varies slightly but this variation has been 
neglected. For the second stage, untrimmed aerodynam- 
ics is used so the eg position is not needed. 

Propulsion 

The ALS is powered by ten liquid hydrogen/liquid 
oxygen low cost rocket engines (LCE): seven power the 
booster and three power the core. All engines are ig- 
nited at launch; staging occurs when the booster fuel is 
depleted; and the core engines burn until insertion. 

Propulsion characteristics of interest are the thrust 
T, vacuum thrust T vac > and the specific impulse J ip (see 
Eqs. 1). If the exit pressure is conservatively approxi- 
mated as p e = 0, the thrust of a single engine is modeled 
as 

T* = T vac 9 — pA t 9 (12) 

where the prime denotes one engine, p is the atmospheric 
pressure at the altitude of the rocket, and A% is the exit 
area. Date relevant to one LCE are as follows: 

T vac ' = 580, 110.0 lb 
A t * s 40.381 ft 2 (13) 

V = 430.0 sec . 


Vm =Cm(M,o) U«) 


where M denotes the Mach number and the bar indi- 
cates that the moment is about a fixed point (launch 
eg). About the actual center of gravity, the moment is 
given by 


Cm — Cm Cj) 


Vcg - 10.36 
/ 


(17) 


since x C8 is assumed not to change. 

While the aerodynamic data could have been used in 
tabular form with linear interpolation to read the tables, 
the approach taken is to assume polynomials in a with 
Mach-number-dependent coefficients. For the first stage, 
the coefficients are written as 


C D = Cdo(M) + C DoJ (M)a s + C D „ s (M)a 3 
C L = ( 18 ) 

73 m — 73 mo (M) + 5m, (M)ct 


where the Mach-number-dependent terms have been ob- 
tained from cubic-spline curve fits of the tabular data. 
After staging, the flow regime is hypersonic and the aero- 
dynamic force coefficients are modeled as 

Cj) = Cd 0 + Cn^ct + Cd^ol 1 

Ci = Ci m <x + Ci m2 or 3 (10) 


For the complete vehicle, 

T — kT* , I$p = I$p i — kT V ae (1^) 

where k = 10 before staging and k = 3 after staging. 
Specific impulse is like specific propellant consumption 
(weight flow rate of propellant per pound of thrust); 
hence, it has the same value regardless of the number 
of engines operating. 


where the coefficients of a are constants. Also, pitching 
moments are assumed to be negligible after staging, that 
is, untrimmed aerodynamics are used (6 = 0). 

A peculiarity of the aerodynamics of the combined 
vehicle at supersonic and hypersonic speeds is that the 
drag coefficient has a minimum at a positive angle of 
attack. This is caused by the aerodynamic shading of 
the booster by the flow field of the core. 


Aerodynamics 

The drag, lift, and pitching moment are related to 
their respective coefficients by the standard equations 

D = qSbCv , L ~ qSbCL = qS^C m (15) 

where q = \pV 2 is the dynamic pressure, St = 
1413.71 ft 2 is the cross-sectional area of the combined 
vehicle (booster 4- core), and / is the length of the core. 
While the aerodynamic coefficients are needed at and 
about the center of gravity (eg), the aerodynamic data 
has been provided at and about the launch eg. Although 
the drag and lift transfer directly, the moment changes 
with eg position. Therefore, the aerodynamic data at 
the eg must be related to the launch eg. 

The aerodynamic data are preliminary estimates as- 
sociated with the development of the six-degreee-of- 
freedom simulation presented in Ref. 4. These data are 
provided in tabular form (Tables 2 through 6) consistent 
with the functional relations 

C D = C D {M,a) = Ci(My<x) , 


III. The Optimal Control Problem 

Formally the optimal control problem considered here 
is to find the control history u(t) which minimizes a per- 
formance index of the form 

J = <*(*;) (20) 

subject to the differential constraints 

z = /(x,u), (21) 

the prescribed initial conditions 

to = to. » x o = %o, » (22) 

the prescribed final conditions 

*(*/) = 0 , (23) 

and a state-variable inequality constraint 

S{x) < 0 . (24) 
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Each of these quantities is discussed below. 


State Variables and Control Varables 

The state variables are x T = [A r h V y tp m] while 
the control variables are u T = [a /*]. 

Performance Index 


It is desired to maximize the final mass. Hence, the 
performance index is taken to be 


$ = - 



(25) 


where the minus sign is included because the perfor- 
mance index is actually minimized and where m rt j is 
the sum of the payload mass, the payload margin mass, 
and the payload fairing mass. A performance index of 
$ = —1.0 means that the reference mass is inserted into 
orbit with no extra fuel. 


Differential Constraints 


The differential constraints are the equations of mo- 
tion (Eqs. 1) completely expressed in terms of the state 
variables and the control variables. 

Prescribed Initial Conditions 

For the trajectory design problem, the initial condi- 
tions are taken to be 

t 0 = Osec , A 0 = -80.54 deg ,r 0 = 28.5 deg 

h 0 = Oft , V„ = 0 — ,70 — 90.0 deg (26) 
sec 

tp 0 = 0.0 deg ,m 0 = 108, 574.93 slugs 

During the vertical rise segment, the heading angle is 
undefined, so the initial condition on tp is actually the 
heading angle during the pitch-over segment. 


Prescribed Final Conditions 


The Advanced Launch System is being designed to 
place a nominal payload at perigee of an 80nm by 150nm 
transfer orbit of 28.5 deg inclination. As a consequence, 
the equality constraint residuals are 

= h, -486,080ft , = V it -25,776.9 — 

J sec 

= 7/ , ¥4 = cos if — cos(28.5 deg) (27) 

where the inertial velocity and the inclination are related 
to the relative states as follows (Ref. 5 and 6) 

Vi = [ V 2 + 2Vrwco$ycosxpcosr -f (rwcosr) 2 ] ^ (28) 


cost = 


cost (V cosy cos + rwcosr) 

[V 2 cos 2 y + 2Vru> cosy cosx^co sr -f (rucosr) 2 ]i 


State- Variable Inequality Constraint 

Based on structural considerations, the ALS must 
not exceed a maximum dynamic pressure of q = 
650 lb/ft 2 . Therefore, the state- variable inequality con- 
straint residue S is 

s = \pV 2 - 650 lb/ft 2 . (29) 

* 

Actually, in a standard atmosphere, the limit is q ma s = 
850 lb/ft 2 . The value of 650 lb/ft 2 is chosen because the 
value of p is approximately 20% smaller in the exponen- 
tial atmosphere than the standard atmosphere around 
the maximum dynamic pressure portion of the trajec- 
tory. 


TV, The Suboptimal Control Problem 


The optimal control problem is converted to a param- 
eter optimization problem (suboptimal control problem) 
as follow: (a) the time is normalized by introducing the 
transformation r = ~\ (b) the control ti(<) is replaced 
by a set of nodal points which is linearly interpolated, 
and (c) the state-variable inequality constraint is con- 
verted to a point constraint by using a penalty function. 

Because of the time transformation, the boundary val- 
ues of r are given by 

_ 3 4 

T 0 — 0 , Tp — 7“ j T\ — , 

*/ */ 


153.54 

r, = — , r f = 1 


(30) 


where t p = 3 sec is the time at the beginning of pitch- 
over and <1 = 4 sec is the time when three-dimensional 
flight begins. Staging occurs when all of the booster 
propellant is consumed; hence, t s = 153.54 sec. 

Figure 2 shows the arrangement of nodal points in 
each stage. Nine nodes are used for the control during 
the first stage, and five for the control during the second 
stage. Even though the duration of the first stage is 
shorter than that of the second, there is more activity in 
a during the first stage, making more nodes desireable. 
The nodes are equally sp diced in each stage so that the 
node times are 


n 


T\ + 



(»•- i) . 1'= 1 
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tv = r,+ I— ^(« - 10) , i = 10 - 14 . (31 ) 
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u (rad) 


o 



Figure 2: Example Control History 


MU- 


(38) 

= 0 

and the inequality constraint 

C 6 = Pj{X) > 0 . (39) 

Derivatives required by the nonlinear programming al- 
gorithm are computed by central differences. 


^ v- 

v 'f. 

c 3 = 7/(*)=0 

_ cosif(X) _ 1 


Note that there are two control nodes at the stage time. 
This has been done in order to find the true suboptimal 
control. 

The dynamic pressure constraint is converted to a pa- 
rameter inequality constraint by introducing the penalty 
function 

P = - f* m*n 2 [(l - -i-),0 ]dt > 0 (32) 

Ju 9majf 

which accumulates value when q > 9 m «** The constraint 
becomes 

Pj> 0 . (33) 

To compute Pj, the penalty function is differentiated to 

form a 

P = — mtn 2 [(l - — ), 0] (34) 

9mox 

where 

P 0 = 0 . (35) 

In all, the nonlinear programming problem involves 30 
parameters, that is, the parameter vector is given by 

X = [fl , on .... , <*u . Mi .••• » Mn . */ ] (36) 


where 6 is the pitch rate during the pitch-over, a* , /it 
are the angle of attack and the bank angle nodes, and 
l] is the final time. 

If values of the parameters (36) are known, the differ- 
ential equations (1) and (34) can be integrated through 
the mission to determine the states and P at the fi- 
nal time. Then, the performance index (25), the orbital 
insertion equality constraint residuals (27), and the dy- 
namic pressure inequality constraint (29) can be com- 
puted. It follows that the performance index and the 
constraints are functions of the parameters (36) such 
that the nonlinear programming problem can be ex- 
pressed as follows: 

Find the set of parameters X which minimizes the per- 
formance index 

j _ _ m /.(£l (37) 

m rt ) 

subject to the equality constraints 

c, . M2 -i = o 

h l. 


V. Numerical Results 

The optimal trajectory has been computed using a 
nonlinear programming code known as VF02AD which is 
based on quadratic programming. Optimal control his- 
tories are presented in Fig. 3, while the resulting states 
are shown in Figs. 4 through 7. The magnitude of the 
performance index is 103.94% where 100% = 171,120 
lb. This means that an additional 6,742 lb of payload 
can be placed in orbit with this vehicle by using the op- 
timal trajectory. The vehicle is inserted into orbit at 
tj = 363.8 sec and the optimal value of the pitch rate 
during the 1.0 sec pitch-over is -.02005 rad/sec. 

Shown in Fig. 8 is the dynamic pressure. It is seen 
that the maximum dynamic pressure occurs at a single 
point and not along a q max subarc. This is due to the 
no-throttling design of the vehicle and the fact that the 
aerodynamic forces needed to fly along q = q m ox cannot 
be achieved. Optimal trajectories with lower values of 
q mag have been calculated, and the results are the same. 

It is difficult to completely determine the meanings of 
the optimal control histories because performance-index 
minimization and constraint satisfaction are going on 
all through the trajectory. For angle of attack, it is seen 
from Fig. 3 that the vehicle initially goes to positive a to 
achieve altitude and decrease q. Then, the dip in a from 
t = 40 to 60 sec allows the vehicle to pass through the 
transonic regime efficiently (Mach 1 occurs at t « 50 sec) 
and to satisfy the dynamic pressure inequality constraint 
(qmax occurs at t « 70 sec). Next, the vehicle returns to 
positive a to get low drag and to decrease the magnitude 
of j. Staging occurs around Mach 8 and the roll off in a 
from positive to negative values during the second stage 
helps pull the trajectory down to meet the final condi- 
tions. For the velocity roll angle, the nonzero values at 
the beginning of the trajectory seem to be caused by the 
rotational effects of earth where the vehicle wants to fly 
at constant latitude throughout most of the first stage. 
Changes in /i near the end of the trajectory help cause 
constraint satisfaction, particularly in the orbit inclina- 
tion. j 

Additional optimal trajectories have been computed 
with the intent of determining what kinds of approxi- 
mations can be made in order to obtain approximate 
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analytical solutions for guidance purposes. First, the ef- 
fect of using untrimmed aerodynamics (6 = 0) rather 
than trimmed aerodynamics is shown in Fig. 9 and 10 to 
change only slightly the optimal controls and to cause a 
relative change in the performance index of 0.2% (376.5 
lb). Hence, untrimmed aerodynamics is a reasonable 
approximation. Second, the question of whether or not 
atmospheric effects can be considered a perturbation is 
considered. This means that the pressure term in the 
thrust and the aerodynamics are neglected; however, the 
dynamic pressure constraint is maintained because it is a 
structural constraint. The optimal controls for this case 
are shown in Fig. 1 1 and 12 and lead to a relative increase 
in the performance index of 16% (27,379 lb). Trajectory 
profiles for the atmosphere and no-atmosphere cases are 
shown in Fig. 13. The optimal control which results 
from the no-atmosphere case is reasonably close to that 
of the atmosphere case and has the same general trend. 
This seems to indicate that atmospheric effects can be 
treated as a perturbation. 


there is a difference of 14 deg between the atmosphere 
and no-atmosphere solutions. Since this region consti- 
tutes less than fifteen percent of the whole trajectory, 
treating atmospheric effects as perturbations could yield 
satisfactory results. 
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Figure 3: TVimmed Aerodynamic Control Histories 
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Figure 7: Latitude and Longitude vs. Time; Trimmed 
Flight 
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Figure 5: Flight Path and Heading Angle vs. Time; 

Trimmed Flight Figure 8: Dynamic Pressure vs. Time; Trimmed Flight 
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Figure 6: Altitude and Velocity vs. Time; Trimmed 
Flight 



Figure 9: Angle of Attack vs. Time 
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Figure 10: Velocity Roll Angle vs. Time 



Figure 11: Angle of Attack vs. Time 




Figure 13: Trajectory Profiles; Atmosphere and 
No-atmosphere 


Table 2. Lift Coefficient (core + booster) 


Subsonic Data 




Angle of Attack (deg) 



M 

±0.0 

±2.0 

±4.0 

±6.0 

±8.0 

±10.0 

0.0 

0.0 

0.06876 

0.1775 

0.2663 

0.355 

0.4438 

0.2 

0.0 

0.08876 

0.1775 

0.2663 

0.355 

0.4438 

0.4 

0.0 

0.08876 

0.1775 

0.2663 

0.355 

0.4438 

0.6 

0.0 

0.08876 

0.1775 

0.2663 

0.355 

0.4438 

0.8 

0.0 

0.08876 

0.1775 

0.2663 

0.355 

0.4438 

1.0 

0.0 

0.08720 

0.1744 

0.2616 

0.3488 

0.4360 

Supersonic/Hypersonic Data 






Angle of Attack (deg) 



M 

0.0 

2.0 

4.0 

6.0 

8.0 

10.0 

1.2 

0.0 

0.0862 

0.1724 

0.2586 

0.3448 

0.431 

1.5 

0.0 

0.086 

0.171 

0.260 

0.351 

0.431 

2.0 

0.0 

0.090 

0.175 

0.262 

0.354 

0.435 

2.5 

0.0 

0.098 

0.181 

0.268 

0.370 

0.460 

3.0 

0.0 

0.100 

0.192 

0.278 

0.385 

0.490 

3.5 

0.0 

0.102 

0.200 

0.290 

0.401 

0.510 

4.0 

0.0 

0.104 

0.202 

0.291 

0.405 

0.510 

5.0 

0.0 

0.104 

0.206 

0.298 

0.410 

0.509 

6.0 

0.0 

0.103 

0.203 

0.300 

0.408 

0.508 

7.0 

0.0 

0.100 

0.195 

0.298 

0.400 

0.502 

8.0 

0.0 

0.095 

0.185 

0.290 

0.395 

0.500 


M 

-2.0 

Angle of Attack (deg) 
-4.0 -6.0 -8.0 

-10.0 

1.2 

-0.084 

-0.170 

-0.260 

-0.350 

-0.431 

1.5 

-0.086 

-0.171 

-0.260 

-0.351 

-0.431 

2.0 

-0.090 

-0.175 

-0.262 

-0.354 

-0.435 

2.5 

-0.098 

-0.181 

-0.268 

-0.370 

-0.460 

3.0 

-0.100 

-0.192 

-0.278 

-0.385 

-0.490 

3.5 

-0.120 

-0.200 

-0.290 

-0.401 

-0.510 

4.0 

-0.120 

-0.215 

-0.310 

-0.420 

-0.520 

5.0 

-0.120 

-0.225 

-0.327 

-0.442 

-0.542 

6.0 

-0.125 

-0.225 

-0.334 

-0.451 

-0.567 

7.0 

-0.115 

-0.222 

-0.332 

-0.452 

-0.580 

8.0 

-0.110 

-0.218 

-0.325 

-0.450 

-0.565 


t 


Figure 12: Velocity Roll Angle vs. Time 
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Table 3. Drag Coefficient (core + booster) 


M 


- 2.0 


Angle of Attack (deg) 
- 4.0 - 6.0 


- 8.0 


Subsonic Data 


0.0 

0.2 

0.4 

0.6 

0.8 


- 0.0201 

- 0.020 

-0.019 

-0.018 

-0.016 


-0.044 

-0.044 

-0.043 

-0.042 

-0.040 


-0.067 

-0.067 

-0.067 

-0.066 

-0.064 

-0.051 


-0.091 

-0.091 

-0.091 

-0.089 

-0.087 

-0.075 


- 10.0 

- 0.115 


-0.115 

-0.115 

-0.113 

- 0.111 

-0.098 


M 

±0.0 

±2.0 

±4.0 

±6.0 

±8.0 

±10.0 

1 . V / 
1.2 

0.003 

-0.029 

-0.058 

-0.089 

-0.119 

0.0 

0.1870 

0.1904 

0.2024 

0.2254 

0.262 

0.314 

1.5 

0.009 

-0.019 

-0.048 

-0.077 

-0.106 

0.2 

0.1872 

0.1906 

0.2026 

0.2256 

0.2622 

0.3142 

2.0 

0.009 

-0.0155 

-0.045 

-0.071 

-0.097 

0.4 

0.2062 

0.2096 

0.2216 

0.2446 

0.2812 

0.3340 

2.5 

0.007 

-0.016 

-0.043 

-0.067 

-0.092 

0.6 

0.2599 

0.2633 

0.2753 

0.2983 

0.3349 

0.3877 

3.0 

0.005 

-0.018 

-0.041 

-0.063 

-0.089 

0.8 

0.3480 

0.3514 

0.3634 

0.3864 

0.4230 

0.4758 

3.5 

0.004 

-0.018 

-0.040 

-0.062 

-0.086 

1.0 

0.7800 

0.7834 

0.7954 

0.8184 

0.8550 

0.9078 

4.0 

0.004 

-0.019 

-0.040 

-0.062 

-0.085 


Supersonic/Hypersonic Data 


M 


0.0 


Angle 

2.0 


of Attack (deg) 
4.0 6.0 


8.0 


10.0 


5.0 

6.0 

7.0 

8.0 


0.005 

0.008 

0.008 

0.008 


-0.018 

-0.017 

-0.017 

-0.017 


-0.038 

-0.028 

-0.028 

-0.028 


-0.058 

-0.058 

-0.058 

-0.058 


-0.082 

-0.078 

-0.076 

-0.075 


1.2 0.800 

0.805 

0.815 

0.838 

0.875 

0.928 






1.5 0.740 

0.703 

0.645 

0.640 

0.635 

0.635 


Table 5. 

Lift Coefficient (core 

vehicle) 

2.0 0.672 

0.656 

0.555 

0.525 

0.525 

0.525 






2.5 0,648 

0.628 

0.512 

0.468 

0.465 

0.455 



Angle of Attack (deg) 



3.0 0.637 

0.608 

0.486 

0.448 

0.431 

0.418 

M 

±0.0 

±2.0 ±4.0 ±6.0 

±8.0 

±10.0 

3.5 0.630 

0.596 

0.470 

0.425 

0.406 

0.392 

8.0 

0.2062 

0.2089 0.2206 0.2417 

0.2733 

0.3160 

4.0 0.628 

0.587 

0.460 

0.410 

0.385 

0.368 

10.0 

0.2X80 

0.2201 0.2313 0.2523 

0.2835 

0.3262 

5.0 0.620 

0.572 

0.448 

0.392 

0.355 

0.352 

12.0 

0.2353 

0.2374 0.2486 0.2703 

0.3024 

0.3459 

6.0 0.617 

0.570 

0.446 

0.382 

0,348 

0.348 






7.0 0.615 

0.567 

0.445 

0.378 

0.340 

0.340 

Table 6. 

Drait Coefficient (core vehicle) 

8.0 0.615 

0.565 

0.445 

0.372 

0.340 

0.338 














Angle of Attack (deg) 



Angle of Attack (deg) 



M 

±0.0 

±2.0 ±4.0 ±6.0 

±8.0 

±10.0 

M -2.0 

-4.0 

-6.0 

-8.0 

-10.0 


8.0 

0.2062 

0.2089 0.2206 0.2417 

0.2733 

0.3160 

1.2 0.803 

0.815 

0.838 

0.875 

0.928 


10.0 

0.2180 

0.2201 0.2313 0.2523 

0.2835 

0.3262 

1.5 0.745 

0.750 

0.773 

0.800 

0.871 


12.0 

0.2353 

0.2374 0.2486 0.2703 

0.3024 

0.3459 

2.0 0.690 

0.708 

0.731 

0.768 

0.822 







2.5 0.665 

0,680 

0.706 

0.745 

0.790 







3.0 0.648 

0.651 

0.688 

0.730 

0.771 







3.5 0.640 

0.650 

0.675 

0.716 

0.757 







4.0 0.631 

0.641 

0.665 

0.706 

0.745 







5.0 0.625 

0.635 

0.651 

0.692 

0.731 







6.0 0.610 

0.625 

0.646 

0.686 

0.727 







7.0 0.610 

0.620 

0.640 

0.685 

0.730 







8.0 0.610 

0.620 

0.640 

0.684 

0.725 








Table 4. Pitching Moment Coefficient (core + 
booster) 


Angle of Attack (deg) 


M 

0.0 

2.0 

4.0 

6.0 

0.0 

0.0035 

0.0271 

0.0508 

0.0744 

0.2 

0.0035 

0.0271 

0.0508 

0.0744 

0.4 

0.0040 

0.0276 

0.0513 

0.0745 

0.6 

0.0052 

0.0288 

0.0538 

0.0757 

0.8 

0.0072 

0.0308 

0.0558 

0.0777 

1.0 

0.020 

0.046 

0.072 

0.098 

1.2 

0.033 

0.062 

0.093 

0.123 

1.5 

0.038 

0.066 

0.095 

0.124 

2.0 

0.033 

0.059 

0.085 

0.111 

2.5 

0.030 

0.056 

0.077 

0.097 

3.0 

0.029 

0.052 

0.071 

0.087 

3.5 

0.028 

0.049 

0.066 

0.080 

4.0 

0.027 

0.047 

0.061 

0.076 

5.0 

0.026 

0.045 

0.057 

0.068 

6.0 

0.026 

0.042 

0.054 

0.068 

7.0 

0.0255 

0.042 

0.053 

0.068 

8.0 

0.0255 

0.042 

0.052 

0.068 


8.0 10.0 


0.0981 

0.1217 

0,0981 

0.1217 

0.0986 

0.1222 

0.0998 

0.1234 

0.1018 

0.1254 

0.124 

0.150 

0.153 

0.183 

0.153 

0.182 

0.134 

0.162 

0.120 

0.135 

0.103 

0.116 

0.094 

0.099 

0.0895 

0.099 

0.085 

0.098 

0.082 

0.096 

0.082 

0.096 

0.082 

0.097 


10 



A SHOOTING APPROACH TO SUBOPTIMAL CONTROL PP^ 

David Hull' ami JyhJotig Shorn* 

Department of Aerospace F.nginooring an<l I'nRineering Mcrliamcs 
The University of Texas at Austin 


Abstract 

The shooting method is used to solve the suboptimal control prob- 
lem where the control history is assumed to be piecewise linear. 
Suboptimal solutions can be obtained without difficulty and can 
by increasing the number of nodes lead to accurate approximate 
controls and good starting multipliers for the regular shooting 
method. Optimal planar launch trajectories are presented for the 
Advanced Launch System. 


Formally, this fixed-final-time suboptimal control problem is 

to find the parameters </ which minimize the perfor- 

mance index 

J m Hxf,ti) + £ I(r,z, «!,..• ,Um,tj)dr (4) 

subject to the dynamics 

*'«■ /(r.x.ui * ( 5 ) 


1. Introduction 

The original motivation for using the shooting method to solve the 
suboptimal control problem (piecewise linear control) has been 
to calculate an accurate suboptimal control and ultimately to 
find the corresponding neighboring extremal feedback control rule. 
Since aerospace minima are usually quite flat, an approximate op- 
timal control can deliver most of the optimal performance. Then, 
the ability to compute the suboptimal control and the neighboring 
extremal without difficulty would be useful. 

In this paper, the shooting method is developed for the subop- 
t jm.l control problem and used to optimize the Advanced Launch 
System trajectory. The usual sensitivity of the solution process 
to the initial guesses disappears completely, and solutions are ob- 
tained without difficulty. Of course, only an approximate optimal 
control is achieved, but if it is not good enought, its accuracy can 
be improved by increasing the number of control nodes. 


2. Suboptimal Control Problem 

The standard optimal control problem is to find the control u(f) 
which minimizes the Scalar performance index 

J m *(*/,</) + /*' 0) 

subject to the system dynamics 


the prescribed boundary conditions 

Tb « 0 , *o * *o.i T J m ** 0 * ( 6 ) 

In these equations, the prime denotes a derivative with respect to 


r, and 


where 


i(r,*,u, «»</) - I 

/(T,*,u lt ... ,«•,*/) - <//(</'.*.«) 

; (r) n<r£n*i 
v ' U+l - Tk 


(7) 


( 8 ) 


and the node times q are fixed. , 

By the usual arguments of the calculus of variations, the equa- 
tions defining the suboptimal solution are given by 


(9) 


and 


xfmj 

Xm-BJ 

RmL + \ T L 

f’R M d r -0 

km l,...,ro 

tit, dr * -Gti 



*o * r / * 1 » 

4 > m 0 , A/ a 


( 10 ) 


X *=/(<, x,u) , (2) 

and the prescribed boundary conditions 

io« 0 , *o “ *o. . ( 3 ) 

The dimensions of x, u, and tj> are n x 1, r x 1, and p x 1, respec- 
tively. This problem is made into a suboptimal control problem 
by normalizing the final time through the transformation r - </«/ 
and by restricting the class of functions to which the optimal con- 
trol can belong. Here, the restricted class is that of piecewise 
linear functions The end points U|, .... «• of the straight line 
segments are called nodes. 


Shooting Mithod 

, ut the suboptimal control problem In a form suitable for ap- 
j, the shooting method, new states v„(r) and u>(r) are m- 
uced to eliminate the integrals in Eq. (9). The optimality 


and 


«*-/. 

A'.-j/J 



km l. 


,m 


(U) 


*o =• *o.. v*. * £, W» ■ 0 

=* 0, A/ = G> ; , 

v k , =0, w, = -G, t . 


1 M . J . Thompaon RegenU ProfcMor 
’Graduate Jteaeareh Aasiatant 


r o »0, 

T/=*i, 


( 12 ) 



5. Conclusions 


If a nrw staU’ vrclor c(0 iMbu'tl a- s 

S = [r T \ T v % ... r m n-\ (”) 

an. I a parameter vector is introduced as 

<i r = (Ul . . . Urn </l . (' 4 ) 

the differential equations (11) can be rewritten in the form 

z = F(r,*,o) ( 15 ) 

Of the initial states, there are 1 + n + m + 1 conditions; only 
X 0 is unknown. At the final time, there are in Eqs. (12) 1 + P + 

„ + m + 1 final conditions. Of these, p equations are solved for 
the p Lagrange multipliers v which are in turn eliminated from 
the remaining conditions to form 

A(z/,o) * 0 (16) 

whose dimension is (n + m + 1) X 1. _ . 

The derivation of the equations for the shooting method is 
straightforward and leads to the following algorithm: 

1. Guess Ao and a 

2. Integrate from r 0 = 0 to T/ = i 

z 9 s F * 0 known 

$' 2 s F x $i $*> = 10/0 0] T (17) 

¥'*F,# + F* #o = 0 

3. Calculate \\h\\. 

4 . Calculate 6 X 0 and 6 a by solving 

[£]--* m 

and using a norm reduction scheme to determine o. 

5. Check for convergence (||A|| < «)• If not, go to 2. 

The advantage of this method is that there is absolutely no 
influence of Ao on *. On the other hand, the sensitivity of the 
shooting method to Ao is replaced by having to accept an approx- 
imate solution. However, by using a reasonable number of nodes, 
it should be possible to obtain Ao’s for which the exact shooting 
method can be converged. 


4. Optimal Planar Trajectory for the ALS 

The Advanced Launch System is a two-stage rocket consisting of 
a core with a side-mounted booster. Staging occurs at the fixed 
time of burnout of the booster. Ref. 1 contains a description of 
the physical model. 

In the optimization problem, the performance index is the final 
mass; the state equations are the equations of motion lor flight u» 
a great circle plane over a nonrotating spherical earth where the 
control is the angle of attack; the initial conditions are all *^‘"“1 
and final conditions are imposed on altitude, velocity, and flight 

path angle. , 

Converged results are presented in Table 1 and Fig. 1 for three 
first and second stage node arrangements. Starting multipliers for 
the 3-2 case are given in Table 1, and the control nodes are taken 
to be a = -.5, 10., 6., 4., -4. deg and i, = 300. sec. Convergence 
required 17 iterations and 241 sec of CPU time on a CDC Cyber 
computer. Also presented are the converged values obtained from 
the standard shooting method. Note that the optimal results are 
approached as the number of nodes is increased. . 

Other than having to derive the multiplier equations, no diffi- 
culties have been encountered during these calculations. 


The shooting approach to suboplin.al control is an effective way 
to obtain approximate optimal trajectories and to obtain starting 
Lagrange multipliers for the regular shooting method. 
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Table 1: Converged Results 



node pattern I 
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— n 
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0.8529 
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371.63 
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-4.125E-4 

~t~ — 

1.0 
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— 
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1.0 

-2.004E-5 
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-2.013E-5 



Figure 1: Angle of Attack Histories 
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NEIGHBORING EXTREMAL GUIDANCE FOR SYSTEMS / jpc J /> # 
WITH A PIECEWISE LINEAR CONTROL , 
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David G. Hull 1 and Clifford E.Helfrich 2 
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ABSTRACT 

The neighboring extremal feedback control law is de- 
veloped for systems with a piecewise linear control for 
the case where the optimal control is obtained by non- 
linear programming techniques. To develop the control 
perturbation for a given deviation from the nominal 
path, the second variation is minimized subject to the 
constraint that the final conditions be satisfied. This 
process leads to a feedback relationship between the 
control perturbation and the measured deviation from 
the nominal state. A simple example, the lunar launch 
problem, is used to demonstrate the validity of the guid- 
ance law. For model errors on the order of 5%, the 
results indicate that 5% errors occur in the final condi- 
tions. 


INTRODUCTION 

In order to develop the neighboring optimal guidance 
law for a dynamical system, it is first necessary to ob- 
tain the optimal control, and this can be a formidable 
task. Currently, most trajectory optimization is ac- 
complished by restricting the class of control functions 
to some subclass, say piecewise linear functions (sub- 
optimal control). Then, the control variables are pa- 
rameters (nodes of piecewise linear function), and the 
suboptimal control is found by applying nonlinear pro- 
gramming methods. Hence, the subject of this paper 
is the development of the neighboring suboptimal feed- 
back control law, assuming that the suboptimal control 
law is available. 

Given the suboptimal control and a perturbation in 
the state at some time, the neighboring suboptimal con- 
trol is found by minimizing the increase is the perfor- 
mance index subject to the constraint that the final 
conditions must be satisfied. Since the first variation 
vanishes, minimizing the increase in the performance 
index is equivalent to minimizing the second variation. 


*M. J. Thompson Regents Professor, Associate Fellow AIAA 
s Graduate Research Assistant 
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The constraint of satisfying the final conditions is ob- 
tained through the use of transition matrices to the final 
point. The above process leads to an analytical expres- 
sion for the gains of the neighboring suboptimal feed- 
back control law. Because of the simplicity of the con- 
trol law, the suboptimal control rule can be applied to 
the vehicles rather than sample and hold. This should 
allow the sample time to be increased, if errors do not 
grow too rapidly. 

To test this guidance rule, it is applied to a simple 
trajectory problem with various levels of modeling er- 
rors. The results indicate that this guidance approach 
has merit. 


SUBOPTIMAL CQ1 


CUV ft* 


L PROBLEM 


The optimal control problem [1] being considered 
here is to find the control history u(t) which minimizes 
the performance index 

J =z (1) 

subject to the state differential equations 

x = f(t,x,u), (2) 


the prescribed initial conditions 

to = » *0 = *0. i (^) 

and the prescribed final conditions 

Here, this problem is converted into a suboptimal con- 
trol problem [2] by assuming that the controls are piece- 
wise linear, meaning that the unknowns become the 
junction points (nodes) of the linear control segments 
and the final time. 

If a denotes the unknown parameter vector, that is, 
a T = [tj, till, « 12 , .... ti 2 i, «22. ... ]. the suboptimal 
control problem is stated as follows: 

Find the set of parameters a which minimizes the 
performance index 

J = F(a) (■*>) 
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subject to the equality constraints 

C(a) = 0 . (6) 

The differential constraints are an integral part of defin- 
ing the functions F and C and are written as 

dx / \ 

— = g(T,x,a) (7) 

To = 0, — ®0, t T f 1 

where r = t/tj and xo. are the specified values of the 
initial states. 

It is assumed that this problem is solved numerically 
by using a nonlinear programming code, and the next 
step is to find the neighboring suboptimal feedback con- 
trol law. 

NF.TOHRORING SUBOPTIMAL CONTROL 

The solution of the suboptimal control problem gives 
nominal control and state histories to be followed by the 
vehicle. However, because of modelling errors, the ve- 
hicle when using the nominal control deviates from the 
nominal state. Hence, it is desired to find the neighbor- 
ing suboptimal control perturbation which enables the 
vehicle to operate in the neighborhood of the nominal 
trajectory. The general philosophy is to find the con- 
trol perturbation which minimizes the increase in the 
performance index while satisfying the prescribed final 
conditions. 

Since the first variation vanishes along the subopti- 
mal path, the increase in the performance index is the 
second variation 


A J = \6a T G aa 6a (8) 

where G = $ 4- v T 'i is the augmented performance in- 
dex and v is a constant Lagrange multiplier. Once the 
suboptimal control has been obtained numerically, the 
second derivative matrix G aa can be computed numer- 
ically. The next step is to find the constraints on 6o 
which guarantee satisfaction of the final conditions (4). 

The variation of the state equation (7) leads to the 
differential equation 


—6x = g x 6x + g a 6a 

dr 


( 9 ) 


which must be solved subject to the boundary condi- 
tions 


To = To, , 
TJ = 1 , 


6x o = 6x o, 

+ i '> T ,6t j = 0 


Next, the solution of Eq. (9) is assumed to have the 
transition matrix form 

6x = $6xj + 96a ( 11 ) 

where 

<*, = /, */ = 0 (12) 

to guarantee that 6xj — 6xj. Then, substituting Eq. 
(11) into Eq. (9) and equating like coefficients leads to 
the following differential equations: 

* = *** (13) 

9 = g x 9 + g a 

which must be solved subject to the boundary condi- 
tions (12). Once $ and ¥ have been obtained, Eq. (11) 
can be used. 

To satisfy the final condition (10), Eq. (11) is rewrit- 


ten as 


Sxj = 

(14) 

Then, assuming 9t f = 0, Eq. (10) leads to 


xl> Xf $~ l 6x — r/ $“ 1 4 r ^a = 0 . 

(15) 

Applied to r 0 , this equation becomes 


tl> Xf *o l *o6a - : ' 6x <> = 0 

(16) 


and is the constraint on the control node perturbation 
6 a imposed by the final condition. 

The last step is to minimize A J as given by Eq. (8) 
with respect to 6a subject to the constraint (16). Stan- 
dard parameter optimization methods lead to 

6a — K 0 6x 0 (17) 


where the gain K 0 is given by 

Ko = G~^J^o T rpl t - 


(18) 


If .the sampling is performed continuously, the param- 
eter perturbation becomes 

6a- Khz (19) 


where 


K 


G~ a 1 9 T < * %• 


( 20 ) 


These gains can be computed at several values of t 
and stored in the onboard computer for interpolation 


( 10 ) 


2 


purposes. 



Two difficulties occur in the use of Eq. (19) as a 
guidance law. First, ♦ goes to zero as r approaches 
unity so that the computation of the gains becomes 
indefinite (zero over zero). This has been handled in 
the following application by computing the gains at r = 
.950 and r = .975 and extrapolating them to r = 1. 
The second problem is determining the value of r on 
the perturbed path since the perturbed final time is 
unknown. This has been accomplished iteratively by 
guessing 6tj, computing r = t/(tj + 6tj), computing 
6a and, hence, 6t j , and repeating the computation until 
the computed 6t; nearly equals the guessed 6tj. 

FVAMPT/E - LUNAR LAUNCH PROBLEM 

The lunar launch problem has been selected as a sim- 
ple example to illustrate the application of this guidance 
law. The optimal control problem is to find the thrust 
inclination history 0(f) which minimizes the time to in- 
sertion 

J = tj (21) 

subject to the differential constraints 


(i 
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Table 1: Results for 5% Modeling Error in Thrust 
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Table 2: Results for 5% Modeling Error in Gravity 


x = ti 

y = v 

U = Or COS 0 
v = asin0 — g , 

the prescribed initial conditions 


t 0 = z 0 = po = u 0 = Vq = 0 , 


( 22 ) 


(23) 


and the prescribed final conditions 

yo = 50, 000 ft, u 0 = 5, 444 ft/sec, v 0 = 0 ft/sec. 

(24) 

The quantities a and g are the constant thrust acceler- 
ation and lunar acceleration of gravity whose nominal 
values are a = 20.8 ft/sec^ and g — 5.32 ft/sec . 

Using five nodes for the suboptima] control calcula- 
tion leads to 

tj ss 272.7 sec, 

8x = 26.09 deg, 

$2 = 20.68 deg (251 

0 3 = 15.34 deg, V ' 

0 4 = 9.061 deg, 

0 5 = 3.113 deg 

To test the guidance law, a 5% error is introduced 
in a which drives the vehicle away from the nominal. 
Gains have been computed stored at each .025 in r. 
Two implementations have been performed: one is to 


use sample and hold and the other is to use the actual 
linear control. Results are shown for a 4 sec sample 
time in Thble 1. Note that a 5% error in a leads to 
roughly a 5% error in the insertion conditions. 

That the linear control does not do uniformly better 
than sample and hold is disappointing. It is felt that 
the sample time could be increased substantially for the 
linear control relative to sample and hold and still yield 
good results. At any rate these are preliminary results 
and further study is warranted. 

Similar results have been developed for a 5% error in 
g and are shown in Table 2. Qualitatively, these results 
are similar to those in Table 1. 


DISCUSSION AND CONCLUSIONS 

The neighboring extremal feedback control law has 
been developed for systems with a piecewise-linear con- 
trol whose nominal control and trajectory have been 
computed using nonlinear programming techniques. 
Given a perturbation in the state, the neighboring ex- 
tremal control perturbation is obtained by minimizing 
the increase in the performance index relative to the 
nominal value subject to the constraint that the final 
conditions be satisfied. Numerical results for the lunar 


3 































launch problem with mismatches in the thrust accel- 
eration and gravity acceleration show that 5% model 
errors lead to 5% final condition errors. Further study 
of this guidance law seems warranted. 
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Abstract 

The neighboring extremal feedback control law is 
developed for systems with a piecewise linear control 
for the case where the optimal control is obtained by 
nonlinear programming techniques. To develop the 
control perturbation for a given deviation from the 
nominal path, the second variation is minimised sub- 
ject to the constraint that the final conditions be sat- 
isfied. This process leads to a feedback relationship 
between the control perturbation and the measured 
deviation from the nominal state. 

Introduction 

In order to develop the neighboring optimal guid- 
ance law for a dynamical system, it is first neces- 
sary to obtain the optimal control. Currently, moot 
trajectory optimisation (see Ref. 1 for example) is 
accomplished by restricting the class of control func- 
tions to some subclass, say piecewise linear functions 
(suboptimal control). Then, the control variables are 
parameters (nodes of piecewise linear function), and 
the suboptimal control is found by applying nonlin- 
ear programming methods. Hence, the subject of this 
paper is the development of the neighboring subop- 
tirrml feedback control law, assuming that the sub- 
optimal control law is available. 

Suboptimal Control Problem 

The optimal control problem being considered here 
is to find the control history ti(r) which minimises 
the performance index 

j = ( 1 ) 

subject to the state differential equations 

^ = /(r, i,u,</) , (2) 

the prescribed initial conditions 

T 0 ~ T 0., *0 = * 0 ., ( 3 ) 


and the prescribed final conditions 

Tj = 1 , = 0 . ( 4 ) 

Here, the time has been normalized by the final time, 
that is, r = t/tj where tj is an unknown parame- 
ter. This optimal control problem is converted into a 
suboptimal control problem (parameter optimisation 
problem) by assuming that controls are piecewise lin- 
ear, meaning that the unknowns become the nodes 
of the linear control segments and the final time. 

If a denotes the unknown parameter vector, that 
is, a T = [f/, un, uij, . . . , «ji, ujj, . . . ], the differ- 
ential equations (2) and its boundary conditions can 
be rewritten as 

dtx __ | 

— - = f(T,x,o), r 0 = r 0 ., *o = *o., T J - 1 • 

dT (5) 

Given a, these equations can be integrated to obtain 
tj = */(«) so that 4 = = F{a) and 

V> = V[c/(a),</] = C(a) Then, the suboptimal con- 
trol problem is to find the parameter vector a which 
miainruae* the performance index J = F(a) subject 
to the constraint C(a) = 0. 

To solve the suboptimal control problem ana- 
lytically, the augmented performance index J' = 
F( a ) + v T C(a) = G(a, v) is formed. The first vari- 
ation conditions are G m = 0 and C — 0 which de- 
termine o and v. The second variation becomes 
6 7 J' = 6a T G„6a > 0 where CM = 0. 6a can be 
divided into dependent and independent parts, and 
the second variation condition becomes the positive 
definiteness of a matrix. 

At this point, it is assumed that the suboptimal 
control problem is solved by using a nonlinear pro- 
gramming code (see Ref. 1, for example), and the 
next step is to find the neighboring suboptimal con- 
trol. 

Neighboring Suboptimal Control 

The solution of the suboptimal control problem 
gives nominal control and state histories to be fol- 
lowed by the vehicle. However, because of modelling 
errors, the vehicle when using the nominal control 
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deviates from the nominal state. Hence, it is desired 
to find the neighboring suboptima] control pertur- 
bation which enables the vehicle to operate in the 
neighborhood of the nominal trajectory. The gen- 
eral philosophy is to find the control perturbation 
which minimizes the increase in the performance in- 
dex while satisfying the prescribed final conditions. 

Since the first variation vanishes along the subop- 
timal path, the increase in the performance index is 
the second variation 

A J = \ 6a T Gaa6d (6) 


and is the constraint on the control node perturba- 
tion 6a imposed by the final condition. 

The last step is to minimize A J as given by Eq. 
(6) with respect to 6a subject to the constraint (14). 
Standard parameter optimization methods lead to 

6a = Kq6xo (^5) 

where the gain K 0 is given by 


subject to C a 6a = 0 which is imposed below. Once 
the suboptimal control has been obtained, the second 
derivative matrix G«« can be computed numerically. 
The next step is to find the constraints on 6a which 
guarantee satisfaction of the final conditions (4). 

The variation of the state equation (5) leads to the 
differential equation 

y-6z = g t 6z + g a 6a (7) 

aT 

which must be solved subject to the boundary con- 
ditions 

TO = To,, 6z 0 = 6zo. /g\ 

T/ = 1 , 1>t,6xj + J>t,6i] = 0 . 

Next, the solution of Eq. (7) is assumed to have the 
transition matrix form 

6x = *6zj + *6a (9) 


where 

*/ = /, *,=0 (10) 

to guarantee that 6xj — 6zj. Then, substituting Eq. 
(9) into Eq. (7) and equating like coefficients leads 
to the following differential equations: 


*' = 9s* 

*' = 9s* + 9a 


( 11 ) 


which must be solved subject to the boundary con- 
ditions (10). Once * and * have been obtained, Eq. 
(9) can be used. 

To satisfy the final condition (8), Eq. (9) is rewrit- 
ten as 

6z f =*~Hx -*~ l *6a (12) 

Then, for the case where il> tj — 0i Eq. (8) leads to 

V> r/ *" 1 «x-V'r / ^~ 1 *^ = 0 • (13) 

Applied to tq, this equation becomes 

V> r/ *o 1 <Ma - ipx t *o l6x <> = 0 ( 14 ) 


Application 

In Ref. 2, neighboring suboptimal control has been 
applied in the same manner as neighboring optimal 
control, that is, sampling is assumed to occur contin- 
uously so that To = r. However, in optimal control, 
any part of an optimal trajectory to the final con- 
straint manifold is an optimal trajectory, but this is 
not the case in suboptimal control. In fact, there 
may not even be enough nodes between the sample 
p oi nt and the final constraint manifold to satisfy the 
boundary conditions. 

Two alternate approaches are being considered. 
First, additional nodes are placed near the final 
constraint manifold to make neighboring suboptimal 
control valid near the end of the trajectory. Second, 
the suboptimal control is computed from each node 
to the final constraint manifold, and the gains (16) 
are computed at each node. These gains are linearly 
interpolated for the operation of the vehicle. Unfor- 
tunately, no results for either case are available at 
the time of this writing. 
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